Towards an anisotropic atom-atom model for the crystalline 
phases of the molecular Ss compound 

C. Pastorino and Z. Gamba 
Department of Physics, Comision Nacional de Energia Atomica-CAC, 
Av. Libertador 8250, (1429) Buenos Aires, Argentina^ 

Abstract 

We analize two anisotropic atom-atom models used to describe the crystalline a, (3 and 7 phases 
of Ss crystals, the most stable compound of elemental sulfur in solid phases, at ambient pressure 
and T<400 K. The calculations are performed via a series of classical molecular dynamics (MD) 
simulations, with flexible molecular models and using a constant pressure-constant temperature 
algorithm for the numerical simulations. All intramolecular modes that mix with lattice modes, 
and are therefore relevant on the onset of structural phase transitions, are taken into account. 
Comparisons with experimental data and previous results obtained with an isotropic atom-atom 
molecular model are also performed. 
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I. INTRODUCTION 



Our knowledge of the complex phase diagram of elemental sulfur, as well as of the 
statical, dynamical and thermodynamical properties of all its condensed phases, is far from 
complete [|I|, 0]. Its complexity is due to a rich variety of molecular allotropes (open and 
cyclic chains of different lengths) that correlate with a large variety of chemical and physical 
properties in the crystalline and liquid phases |1|, [2j . 

The most stable, and the most studied, sulfur allotrope at ambient temperature and 
pressure (STP) is S 8 , a crown-shaped cyclic molecule. Three crystalline phases, a-, (3- and 
7-Sg, have been reported at ambient pressure and T<400 K. Figs, la, b and c show their 
structures. The orthorhombic a-S 8 can be measured || |J at temperatures up to 385.8 
K, its metastable melting point. At 369 K there is a solid-solid a-f3 phase transition, but 
it is most often observed when the samples are studied from high to low temperatures. 
Monoclinic [3-Ss f| is the most stable phase in the 369-393 K temperature range. /3-S 8 
crystals are usually obtained from the melt Jl[ and they show an orientational order- disorder 
phase transition at 198 K |. A third crystalline allotrope 7-S 8 ]7) (Fig. lc) has been 
observed, with a density, at STP, 5.8 % higher than that of a-S 8 . 



Most theoretical studies of the S 8 crystalline phases, up to very recently, have been 
limited to calculations of the intramolecular mode frequencies for the isolated S 8 molecule 
and their splitting when they are embedded in a a-S 8 matrix, including calculations of their 
infrared and Raman intensities ||. Constant- volume molecular dynamics (MD) simulations 
with the same aim have been also performed f9j. 



We recently reported jL0[ [TT[| a series of constant pressure-constant temperature classical 
MD simulations of a—, (3— and 7 — S% crystals using an extremely simple and flexible 
molecular model, a slight modification of the one proposed in the constant volume MD 
simulations of a-S 8 ||. The simple intermolecular potential was of the Lennard- Jones (LJ) 
atom-atom type, except that its parameters were fitted in ref. |1(] to the experimental 
data on the STP crystalline structure of a-S 8 and an estimated value of its configurational 
energy, calculated as in refs. [H], [11] and based on the experimental STP value of the heat 
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of sublimation of a-Ss |12| (it is also explained in the Calculations section of this paper). 
These LJ parameters (Ei SO = 1.70 kJ/mol and <Xj SO =3.39 A) were afterwards used throughout 
the calculations of all samples of a—, (3— and 7 — S$ crystals at different temperatures 
and zero kbar pressure |TI], |H| . The intramolecular potential included a harmonic well 
term for the SSS angle and a torsional potential of the double well type for the sulfur 
chain, different from hydrocarbon chains that show a triple well torsional potential. These 
intramolecular coordinates fully describe all the low frequency intramolecular modes that 
mix with lattice modes || [J and can, therefore, be relevant for the lattice stability or in 
the onset of structural phase transitions. 



Our calculations showed that the simple and flexible molecule model succesfully repro- 



duced most experimental data available on the three crystals ||10| , pT| , including structural 
and dynamical data of a- and 7-S 8 at STP, the high temperature orient at ionally disordered 
phase of /3-Sg at 375K and its low temperature ordered phase. Nevertheless, below 200 K 
this simple isotropic LJ model did not reproduce the 100 K structure of a-Ss, as determined 
by a neutron and X-ray measurements [13|]. Although no structural phase transition was 
measured for a-Sg crystals ||13|| , when decreasing the temperature of the sample below 200 
K, our orthorhombic sample turned out to be unstable using this isotropic atom-atom 
potential, and transformed to a monoclinic structure that we labeled a'-S 8 |TUL [TTj] (Fig. 1). 



Assuming that the experimental data of ref. [13] do not correspond to a metastable 
state, the problem can be correlated to the extremely simple model molecule and in this 
work we propose to take into account the anisotropy of the sulfur atom due to its 3s 2 
3p 4 electronic configuration. This idea is not new, it gave excellent results in similar 
problems. For example, in the case of alkane hydrocarbon chains it is known that an 
unique isotropic atom representing the CH 2 units is able to reproduce many facts of the 



high temperature phases 14], but agreement with the measured crystalline structures is 



attained only when hydrogen atoms are explicitly included |T5|]. An intermediate model 
was proposed in ref. [16], consisting in displacing an unique isotropic interaction site from 
the C position to the center of mass of the CH 2 unit, this improvement was enough to study 
the liquid phase of alkanes but not their crystalline phases [T7]]. The "all atoms" model 
has afterwards been succesfully applied to many studies, like the conformations of flexible 
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molecules |Tj§ , the equation of state of alkanes |17| and the condensed phases of glycerol |L9 . 



Also in ref. |2(] the parameters of an aspherical atom-atom model were fitted to reproduce 
11 crystalline structures of elemental sulfur, the location of the sites that simulate the 
lone-pair electronic distribution was refined. The calculations were performed in the 
harmonic approximation, for rigid molecules, and no attempt to include temperature, 
intramolecular soft modes, lattice stability or vibrational frequencies was performed. 

In the following sections we present the models of the anisotropic intermolecular potential 
that are analized in the present calculations, the method of calculation, the obtained results 
for a-, a'-, (3— and 7 — 58 crystals and our conclusions. 



II. THE ANISOTROPIC ATOM- ATOM POTENTIAL 

For the present calculations we analize two simple anisotropic atom-atom LJ intermolec- 
ular potential models, based on the experimental charge distribution of a-S$ , where the 
location of the sulfur lone pair orbitals were approximately determined. Our first anisotropic 



case molecule model model of refs. [10], [LI] replaces each isotropic LJ atomic site by four sites 
at a distance d an i from the S atom, two of them located along the SS bonds and the other two 
are generated with the S4 element of symmetry located at the atomic S site. This multi-site 
model (that we called ANI-S4) is a simple way to represent the atomic anisotropy. 
For the distance d an i we considered two values. The first one is d an i=0.Q A, the experimental 



value determined at 100K |13|]. The second is a reduced value d ani =0.1 A, that was included 
after taking into account that our four sites anisotropic model is extremelly simple and do 
not include a site at the S atom, making artificially large the anisotropy of the S atom. 

The new LJ parameters of these sites were adjusted to fit the experimental structure and 
the estimated value of the potential energy of a-S 8 at 300K, obtaining (e ani =0.225 kJ/mol, 
a ani =2.87 A) for d ani =0.Q A, and (e an j=0.115kJ/mol, a an j=3.35 A) for d ani =0.1 A. It has 
to be taken into account that, should we have taken d an i=0A, the LJ parameter's values 
should have been: e an i= £j SO /16=0.107kJ/mol and cr an j=cri SO =3.39 A. 
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As we show in the following sections, due to the results obtained with the ANI-S4 
potential model, for both d an i distances, a further an extremely anisotropic atom model was 
tested, for the sake of completeness. The last model locates the two external sites almost 
perpendicular to the plane determined by the corresponding S atom and the two nearest S 
neighbors, at d an i=0.Q A from the S atom (ANIj_ model). This approximated location is also 
suggested by the experimental data of ref. [TB| and the MO calculations that they include, 
on the sulfur containing simple molecule H2S2. As we explain in the next section, the most 
stable numerical simulations were obtained when the two sites are generated by simplely 
displacing the two sites at a distance d an i=0.6 A from the location of the corresponding 
atom Sj (j=l,8) and orienting them along the vectors l[(S J+2 - S_j+i)+(S_j_2-Sj_i)]. That 
is, the two external sites are located in a plane parallel to that formed by the next two 
nearest SS bonds. This direction is approximately perpendicular to the plane formed by the 
Sj_iSjSj + i atoms and this model is considered a case study that represents an extremely 
large atomic anisotropy. The LJ parameters of the ANIj_ model were also adjusted to fit 
the experimental structure and estimated potential energy of a-S$ at STP, obtaining in this 
case e an i=0.15 kJ/mol and a a ni=3. 17 A. 

The LJ parameters of each multisite model, all of them adjusted to fit the experimental 
STP data of a-Sg, were afterwards used to perform the series of MD runs of the above 
mentioned four crystalline phases, at zero pressure and in the 50K - 450K temperature range. 



The intramolecular terms of the potential model, for all molecular models, are essentially 



equal to those of refs. [JTOJ, |TTJ, with the only difference that the force constants take now a 
value about 10 % lower, in order to improve the fit of the intramolecular frequencies to the 
experimental data. The intramolecular potential for the S-S-S bending angles consists of a 
harmonic well 

V^) = l -C^(3-f3 )\ 

with a force constant of 0^/^=22000 K/rad 2 and (3q= 108 deg., kg is the Boltzmann 
constant. The intramolecular potential for torsion (r) angles is a double well: 

V(t) = A t + B t cos(r) + C T cos 2 (r) + D T cos 3 (r), 
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with Ar/k B =51.473 K, 5 r /k fl =664.573 K, C T /k B =2068.092 K and £> T /k B =476.453 K. 
These parameters describe a double well with minima at r=l 98.8 deg., and a barrier height 
of about 8 kJ/mol at r=180 deg. The barrier height at r=0 deg. is of 27 kJ/mol, out of the 
range of energies explored in these simulations. 

The bond lengths S-S are held constant at their mean experimental value of 2.06 A. 

III. CALCULATIONS 

The method of calculation, as well as the performed series of numerical simulations are 
equal to those included in refs. [H], [□], the only change is in the potential model used. 



Briefly, the phase diagram and dynamical bulk properties of Ss crystals, as given by these 
simple molecular models, are calculated at kbar pressure in the (N,P,T) ensemble, by a 
series of classical MD simulations. The external pressure is considered isotropic and the MD 
algorithm allows volume and shape fluctuations of the sample, in order to balance the inter- 
nal stresses with an external isotropic pressure ETJ . The algorithm controls pressure via an 



extended system which includes as extra variables the MD box parameters, periodic bound- 
ary conditions are applied. The temperature control of the sample follows the approach of 



Nose |22|, |23[, which also includes an external variable. The equations of motion of these 
flexible molecules are integrated using the Verlet algorithm for the atomic displacements and 
the Shake algorithm for the constant bond length constraints on each molecule ||24j| . The MD 



algorithm is based in that used in a study of black Newton films [25], afterwards extended 



to constant-pressure constant-temperature simulations via an extended Lagrangian Eq, 27 



The final version is similar to that used in ref. |28|. In the present work, as the interaction 
sites of our anisotropic atom model are not coincident with the S atomic sites, the algorithm 
employed to translate the forces from massless to massive sites is based on that written for 
rigid molecules|29[, extended to the case of flexible molecules. 

The runs in the (N,P,T) ensemble were performed every 25 or 50 K in the range 50-450 
K. At each point of the phase diagram the samples are equilibrated for 20000 to 30000 time 
steps (each one of O.Olps.) and measured in the following 20 ps. Near the phase transitions 
the stabilization times were increased several times. 

The exception are the series of runs corresponding to the ANI-S4 model with d an i=0.6A, that 
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were performed with a time step of 0.005ps. Nevertheless, these runs, after equilibration 
of the sample at the desired T and P, and in a free trajectory in the phase space, give a 
large deviation of the total energy of the system (including the aditional coordinates of 
the extended system), implying that the numerical algorithm should be improved or that 
the molecular model do not calculate the sample as a stable one. For this model we only 
include, for the sake of comparison, the calculated structure and configurational energy 
that is measured (after the initial transient) in the equilibration trajectory, an algorithm 
that following ref. ^ simulates the system in the canonical ensemble. 

The samples of a-S 8 and ce'-Sg crystals consisted of 3*3*2 orthorhombic cells (288 
molecules). The sample of /3-Sg crystals consisted of 4*4*4 cells (384 molecules) and that 
of 7-Sg crystals of 5*3*5 cells (300 molecules). 



The free energy of all samples is estimated, at each point of the phase diagram, taking into 
account the contribution of the calculated vibrational modes to the entropy of the sample 
|3~I|| . This simple method is widely used in lattice dynamics calculations: In a first step 
the calculated vibrational density of states VDS(z^) is normalized to the total number of 
degrees of freedom of the MD sample. This normalized function VDS(z/) gives the number 
of oscillators within jA_ Au and in a first approximation, by considering a system of harmonic 
oscillators, the calculation of their contribution to the entropy is straightforward In 
the same way the internal energy and enthalpy of the system can be derived, the enthalpy 
being directly comparable to thermodynamic measurements. 



IV. RESULTS 

Using the anisotropic LJ model ANI-S4, with distances d an i=0.1 and 0.6A, the or- 
thorhombic a-Sg sample does not spontaneously show the structural phase transition to 
the monoclinic a'-Sg, as was calculated with the isotropic LJ model when the temperature 
of the orthorhombic sample was decreased below 200 K HTo| , [TT|] . Nevertheless, if a new 
series of simulations are performed on a sample that initially has the structure of the 
monoclinic a'-Sg, this phase is calculated stable, in the time scale of our simulations 
(about 200 ps.), with a potential energy lower than that of a-Sg. Fig. 2 shows the 
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calculated configurational energy (the time and sample average of the potential interaction 
energy per molecule) and volume per molecule as a function of temperature, for the four 
crystalline samples, and using the ANI-S4 model with d an i=0.1 A. Table I includes our 
results for the ANI-S4 model with d ani =0.1 and 0.6 A at selected points of the phase diagram. 

A dramatic change is observed when the ANIj_ potential model is used. These results 
are included in Fig. 3 and they clearly show a large diference with those of Fig. 2. The 
extremely anisotropic case model ANIj_ gives results that, in general, compare better 
with some of the experimental data than those included in Fig. 2. In particular, Fig. 
3 shows that the low temperature orthorhombic a-Sg phase is now calculated with a 
lower configurational energy (and volume per molecule) than the monoclinic a'-Sg phase 
suggested by the isotropic and ANI-S4 models. Nevertheless, this phase is calculated as 
stable with all tested models and this fact implies that a'-Sg might be a metastable phase 
of elemental sulfur. 

As we have also discussed in ref. [Hj, the uncertainty in our comparison with the scarce 
experimental data on (3- and 7-Sg phases is mainly due to contradictory measurements of 
7-Sg crystals. Density measurements at 300K imply that their density is 5.8% higher than 
that of a-Sg at STP 0. On the other hand the X-ray measurements determine a 
crystalline cell with a larger volume per molecule than that of a-Sg, implying a 2.3% lower 
density. Further measurements of this phase are needed. 

The orient at ionally disordered phase of /3-Sg is reproduced by the isotropic [flOfl , the 
ANI-S4 and ANIj_ anisotropic models. At high temperatures, the molecules at (0,0,0) and 
(1/2,1/2,1/2) are dynamically disordered. The inertial axis perpendicular to the molecular 
plane librates around its equilibrium orientation, whilst the molecule reorientates within 
the plane. We calculate a characteristic reorientational time (within molecular plane) of 9 
ps at 375 K with ANI-S 4 {d ani =0.6A) and of about 50 ps at 375 K with ANI_l. Although 
no sharp transition is found for this crystal when decreasing the temperature of the sample, 
below 150 K both molecules are orient at ionally ordered and related by a C2 symmetry 
axes. No measurements on the molecular reorientational times and lattices parameters at 
diferent temperatures are available for comparison. 
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It has to be pointed out that, when analyzing the relative stability of two phases we 
have estimated the Gibbs free energy of our samples by including the contribution of the 
vibrational modes to the entropy Fig. 4 compares the data of a-, a'-, f3- and 7-Sg. At low 
temperatures the relative stability of the phases is equal to that obtained from the potential 
energy data. At high temperatures the merge of all phases is observed. No crossings (which 
identify phase transitions) are found. 

Table I compares the experimental and calculated (using the isotropic and anisotropic 
models) structural data at the starting point of the phase diagram of our series of simula- 
tions. The calculated average pressure, not included in the Table is kbar in all cases. For 
example, for a-Sg at 300K is P=0.0(5) kbar. 



As in refs. [10], [11], the a- (3 phase transition was not found in our calculations, in accord 
with the experimental fact that the transitions are promoted by the disorder existing in the 
real samples 0, if not, metastable states can be maintained for days. As calculated in refs. 



10, 11 with the isotropic model, also for the anisotropic atom-atom model the solid-liquid 



phase transition is located at temperatures near the experimental value (~400 K), but only 
for a disordered sample of Ss molecules. This sample is analized with a MD algorithm that 
allows changes in the MD box volume, but not in its angles, as a function of pressure and 
temperature. The disorder is generated by initially increasing the temperature of a cubic 
sample to 500 K and afterwards following its evolution by increasing and decreasing the 



temperature in the 350-450 K range. The result is totally similar to Fig. 4 of our ref. |1C 
and it is not included here. This last calculation is only of theoretical interest, and we did 
not intend to simulate the real liquid phase of elemental sulfur. The experimental data show 
that the disorder, in the case of the real solid-liquid phase transitions, is generated because 
the Ss molecules start to dissociate when the temperature of the sample is near the melting 
point and our calculations should be valid very near the transition, when the fraction of 
broken molecules is small. 
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V. CONCLUSIONS 



It has to be pointed out that, as we are dealing with molecular crystals, as long as 
no chemical reactions or electronic excitations are involved, the classical approximation 
is good enough for the calculation of structural and dynamical properties, provided 
that the inter- and intramolecular potential gives a good description of the real inter- 
actions. Moreover, simple model molecules can be very helpul to elucidate the complex 
phase diagram of the different crystalline allotropes of elemental sulfur, some of them 
with a large number of molecules per cell and others with molecular reorient at ional disorder. 

In refs. [L0|, |TT] we verified that, surprisingly, a large amount of crystalline properties can 
be very nearly reproduced by a simple and isotropic atom - atom LJ model. The most dis- 
cordant point being the finding of a low temperature monoclinic sample a'-Sg, that has not 
been observed || but, nevertheless, might be a metastable phase of this elemental compound. 

Assuming that this a'-Sg phase is, probably, an artifact of the simple isotropic inter- 
molecular potential model, in this paper we analized two anisotropic case molecular models. 

Our calculations show that an anisotropic ANI-S4 model with local S4 symmetry 
calculates the same stable structures and dynamical properties than the isotropic model, 
in the time scale of our simulations (about 200 ps.), being the isotropic model more sen- 
sitive to changes in the thermodynamic parametes P and T than the anisotropic ANI-S4 one. 

The results obtained with the ANIj_ case model are, instead, qualitatively very dif- 
ferent from the other models. The four crystalline phases are calculated as polymorphic 
structures, but their relative cohesive energy is changed. With this model molecule 
the a-Sg phase has a lower configurational energy than a'-Sg in all the studied range 
of temperatures, in agreement with available measurements. The calculated density of 
7-S 8 crystals is in accord with X-ray measurements, but not with the density data of ref. |2], [7|. 

In conclusion, because of the amount of experimental data that these two case model 
molecules nearly reproduce, and due to the fact that the a'-Sg allotrope is calculated stable 
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(for about 200ps.) with three model potentials, we propose a thorough measurement of 
the crystalline phases of Sg, including different techniques employed to grow the crystals 
at different temperatures. Additional measurements of the (5- and 7-S 8 crystals are also 
needed, in particular measurements of molecular reorientational times for /9-Sg and lattice 
parameters as a function of temperature and pressure. Accurate measurements of the 
density of 7-Sg are still lacking. 

Detailed calculations of the lone pair orientation in elemental sulfur molecules, including 
open chains and cyclic molecules of different size, are presently being carried out using the 
GAMESS (version 2001) ||32] program. Lastly, a further and simple anisotropic model will 
be developed with which we plan to extend our calculations to the crystalline phases of the 
molecular allotropes and the liquid phase of elemental sulfur. 
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Table I: Experimental and calculated lattice parameters: a, b, c (A) and monoclinic 
angle /3(deg.). iso denotes the calculations performed with the isotropic atom-atom 
model fllO|, pTTf] and ANI those performed with the corresponding anisotropic model. Z is 
the number of molecules in the cell, V (A 3 ) is the volume per molecule. See the text for a 
discussion about the 7-S 8 data. 
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Figure: 

Fig. 1: Crystalline phases of the molecular Sg compound. From left to right: a) 
orthorhombic a-Sg at 300 K, b) monoclinic /?-Sg at 150 K (with the high temperature 
orient at ionally disordered molecules in light colour), c) pseudo-hexagonal 7-Sg at 300 K, d) 
other view of the orthorhombic a-Sg at 300 K, e) calculated monoclinic a'-Sg, at 100 K, 



with the isotropic LJ model [110, [TT 



Fig. 2: Configurational energy V con f (kJ/mol) and volume per molecule (A 3 ) at kbar 
and as a function of temperature T (K), as given by the anisotropic ANI-S4 model, with 
d an i = 0-l A. The lines are a guide to the eyes. 

Fig. 3: Configurational energy V con f (kJ/mol) and volume per molecule (A 3 ) at kbar 
and as a function of temperature T (K), as given by the anisotropic ANIj_ model, with 
d a ni=Q-6 A. The lines are a guide to the eyes. 



Fig. 4: Gibbs free energy (kJ/mol) at kbar and as a function of temperature T (K), as 
given by the anisotropic ANIj_ model, with d an i=0.6 A. The lines are a guide to the eyes. 
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This figure "figl.gif" is available in "gif" format from: 



http://arXiv.org/ps/cond-mat/0007059v2 
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